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A two-component Bose-Einstein Condensate 
(TCBEC) exhibits a wide range of interesting be- 
haviour that has been the subject of much theoretical 
and experimental research. The two components may 
form either miscible or immiscible phases exhibiting 
complex density profiles [l|, @, 0, |5[. For repulsive 
interspecies interactions, the two components may also 
form symmetry broken states 0, [H, @, S 8, 9, 10]. For 
rotating TCBECs, interlocking vortex lattices [11] and 
vortex sheets [12] have been predicted theoretically, 
with vortex lattices being confirmed experimentally 
[l3]. This has led to much interest in the dynamics, 
and collective excitations of these systems [13, EH, fl6j . 
and recently a number of papers have also predicted, 
through thermodynamic arguments, the possibility of 
forming giant vortices [17, 18]. 

This large number of phenomena is due to the numer- 
ous experimental parameters that can be varied. For ex- 
ample, the atom number, masses, interaction strengths, 
trapping frequency, and trap ellipticity can all be var- 
ied for each component separately. The parameter space 
is far too large to be fully investigated through numeri- 
cal simulations. For this reason, we investigate the dy- 
namics of a TCBEC through analytic methods then in- 
vestigate points of interest through numerical simula- 
tions. For single component BECs, considerable the- 
oretical [19, 2Q, J21I, [22I, m, [23, [IBI and experimental 
m, 113, 0, 0, (for a summary see Ref. 0) ef- 
fort has been applied to understanding their dynamical 
properties under rotation. A major result of this work 
was that only considering the thermodynamic stability 
of the BEC [45[ does not correctly predict the onset of 
vortex nucleation. It is instead necessary for the BEC to 
be dynamically unstable (dynamical instability implies 
thermodynamic instability, but not visa versa: see, for 
example, Ref. [sij). 

We study the dynamical instabilities of TCBECs in ro- 
tating traps by deriving static solutions in the rotating 
frame, in the Thomas- Fermi limit. We find that these 
solutions describe quadrupolar dynamics. Through nu- 
merical simulations we show that instabilities in these 
solutions lead to phases that have already been pre- 
dicted thermodynamically, such as interlocking vortex 



lattices. In addition, these solutions predict interspecies- 
interaction-mediated centre-of-mass (COM) instabilities. 
We observe these oscillations numerically and find that 
they can settle down into a stable symmetry broken state. 
This state is different topreviously-studied symmetry- 
broken states 0, i, i, 0, II, i, in that it only occurs 
for attractive interspecies interactions. 

The paper is laid out as follows. In Sec. [J we de- 
rive equations that govern the stable irrotational (vor- 
tex free) motion of a TCBEC in a rotating trap. We 
employ the Thomas-Fermi approximation to derive an- 
alytic results and find that when stable, both compo- 
nents undergo quadrupolar oscillations of different mag- 
nitudes. In Seciniwe derive stability equations for these 
static solutions (critical points) in the rotating frame. 
These instabilities indicate that the TCBEC will begin 
to evolve dynamically. We identify four different insta- 
bilities: (i) catastrophic instability, (ii) ripple instability, 
(iii) COM instability, and (iv) intra-species COM insta- 
bility. Types (i) and (ii) lead to turbulence and vortex 
nucleation. Types (iii) and (iv) lead to COM motion. 
In Sec. mil we detail the results of numerical simulations 
used to investigate the analytic predictions. These sim- 
ulations reveal the dynamics induced by the instabilities 
studied in Sec. [Til They show how turbulence in the 
TCBEC allows for the formation of states with topologi- 
cally distinct phase profiles which eventually settle down 
into either giant vortices, interlocking vortex lattices, or 
vortex sheets. The simulations also show that the COM 
instabilities do not lead to vortex nucleation, but even- 
tually settle down into a stable state (in the rotating 
frame), which, for the case of attractive interspecies in- 
teractions, breaks the 180° rotational symmetry of the 
rotating frame Hamiltonian. 



I. THOMAS-FERMI APPROXIMATION 

A. The hydrodynamical equations 

In practice, a TCBEC can be stirred by introducing 
a rotating anisotropy into the confining potential. If 
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done adiabatically, the TCBEC will settle down into a 
state that oscillates in unison with the rotating trap. 
The equations of motion can then be found by consid- 
ering static solutions in the rotating frame. A TCBEC 
can be described by 2 mean-field wavefunctions (^i and 
^2) whose time evolution is dictated by the coupled two- 
component Gross-Pitaevskii Equation (GPE) [32]. In a 
reference frame rotating with angular velocity 11, these 
coupled equations become 
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where the j subscripts (j = 1 or 2) refer to the compo- 
nent under consideration, and / 7^ j. ruj and Vj are the 
mass and the potential affecting component j. gj and 
gi2 are the intra and interspecies interaction coefficients 
given by gj = Anoirfi^aj /ruj and gi2 = 2no7r^^ai2[mi + 
^2]/ [^1^2]. Here, aj and ai2 are the intra and inter- 
species scattering lengths respectively. The no term al- 
lows for a rescaling such that || ||= Nj/riQ where Nj 
is the number of atoms in component j. 

As in the one-component case [19], it is possible to de- 
rive exact solutions in the Thomas- Fermi Approximation 
(TEA) (for a detailed description of the analytic meth- 
ods, see [25]). The GPE in the frame rotating with the 
potential is transformed using = ^/pj(r~t)e'^^^^^''^^ {pj 
is the density of component j, and is the phase), and 
the TEA [33] is applied, giving the hydrodynamical equa- 
tions of motion: 
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1^(0,0,1) and Vj 

+ ujI-z'^]. Here, 
and coj are respectively the ellipticity in the 
x-y plane, the frequency in the z direction, and the 
frequency in the x-y plane when ej = 0. This gives 
a reference frame that is rotating around the z-axis 
in which the potential is static, which corresponds to 
modeling a potential that is rotating around the z-axis 
with angular velocity —Cl. In what follows we enumerate 
the species so that mlujf/gi < m2UJ2/g2' We will see 
later that the behaviour of each component depends 
heavily on this criteria. 

Steady state solutions in the rotating frame are ob- 
tained by setting =0 and —fi^^ = Pj = the chem- 
ical potential of species j. Solving Eq. ([3]) gives two 
possible solutions for the density: if the density of both 
components is non zero, then 
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If one of the components has zero density, then the solu- 
tion for the other component is 
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Regions in the EEC where Eqs. ([5]) and (|6]) are appli- 
cable will be referred to as the overlapping region and 
singular region respectively (indicated by superscripts 
and s). As in the one-component case, the TEA has the 
effect that pj can become negative [33^ . When this is the 
case, we assume that pj =0. Given the above, the total 
density for each component can be expressed as 



p, = pm{p^)H{p'.,) ^ pm{-p-.)H{p]), 



(7) 



where H{p)={) for p< and i^(p)=l for p > 0. Without 
rotation {Vt = V^i = V^2 = 0), Eq. ([7j) corresponds to 
that derived in Refs. [J, Q. 



B. Validity of the Thomas- Fermi Approximation 

In this section we investigate the Thomas-Fermi ap- 
proximation and its validity for a non-rotating TCBEC. 
A two-component BEC has two phases. When —^Jgig2 < 
912 < ^9192 the system is in the miscible phase where 
both components interpenetrate. gi2 > ^gig2 is the im- 
miscible phase [33, '35'] where the two components repel 
each other. In the case of a harmonic trap, the inter- 
species repulsion causes either one component to form a 
shell around the other, or both components to separate 
asymmetrically [4] about the trap centre. The asymmet- 
ric state occurs for gi ^ g2^ ^1 ^ ^2, and gi2 > ^gig2 
[36]. 

For BECs in a harmonic trap, the following general 
behaviour can be seen from Eq. ([5]). When gi2 < 
^2 ^1^1/ (^2(^2) both components form overlapping den- 
sity profiles with concave down parabolic shapes [Fig. [T] 
(I)]. When gi2 = ^2^1^1/(^2^2)5 Pi ^ constant, and, 
for values of gi2 larger than this, begins to dip in 
the overlapping region at the centre of the condensate 
[Fig. [1] (II)]. As gi2 is increased, dips further until 
9i2 > 92P2/P1 where the TEA predicts = at the 
very centre of the trap [Fig. [T] (HI)]. 

The Thomas-Fermi approximation assumes that the 
kinetic energy of the wavefunction is negligible compared 
to the interaction and potential energy terms. For single- 
component BECs, it is valid for large g and atom number, 
which is satisfied for typical experimental parameters. 
These criteria must also be satisfied for each component 
individually if the TEA is applied to a TCBEC. Moreover, 
it has been shown Q that the accuracy of the TEA in the 
two-component case can have a more complicated depen- 
dence on the parameters used. In this section, we explore 
the parameter space of a TCBEC to identify these addi- 
tional constraints. 

We have conducted extensive comparisons between 
the TEA and full GPE solutions, identifying the follow- 
ing conditions for TEA validity: For —^gig2 < gi2 < 
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^2/^2//^i we find excellent agreement [Figs. [T] (I) and (II)] 
{912 = ^2/^2/^1 is the point where the TFA predicts pi{r) 
reaches at r = 0). As gi2 = ^'2^2/^1 is approached the 
TFA results deviate from the full GPE results, but re- 
main qualitatively correct for ^2/^2/^1 < gi2 < \JQ\92 
[Fig. □ (III)]. For gx2 > V^i^ the TFA solutions are 
completely different from the full GPE solutions [Fig. [1] 
(IV)]. This is due to the intense repulsion between the two 
components resulting in sharp curvature in the wavefunc- 
tions, and hence kinetic energy cannot be neglected. For 
912 < —^/gig2 the TFA solutions are unphysical reflect- 
ing the fact that, without kinetic energy, the conden- 
sate is unstable. In what follows we therefore consider 
the TFA solutions only in the region of their validity: 
— ^/9i92 < 912 < \/9i92' In other words, the TFA is 
valid for any two-component BEC that is within the mis- 
cible phase. Eq. ([5]) shows that in this region of validity 
the component with the smallest value of w^uj^ jg (com- 
ponent 1 by definition) will always sit to the outside of 
component 2, and, as expected, this behaviour continues 
into the immiscible phase where component 1 forms an 
almost hollow shell around component 2. 

One should note that the condensates begin to sep- 
arate before the immiscible phase is reached, and the 
transition as g\2 passes ^Jg\g2 is continuous, i.e. there 
is no qualitative difference upon reaching the immisci- 
ble phase. However, as shown by Timmermans [s^, the 
nature of the component separation is different in the 
two phases. For g\2 < \/9\92 any separation is caused 
purely by the potential, if no trap gradient existed then 
both species would overlap with constant densities. For 
9\2 > ^/9i92 the two species separate irrespective of the 
trap. Timmermans points out that this situation is much 
like the case for ordinary fluids. Two fluids may mix 
freely but can still be separated by an external poten- 
tial such as gravity. They will, however, mix freely again 
once stirred. Conversely, immiscible fluids such as oil 
and water always separate, and remain so even after mix- 
ing. Indeed, we find that although there is no observable 
crossover between the two phases while the condensate 
is stable, once the BEC is stirred via rotation, the con- 
densate displays very different behaviour depending on 
which phase it is in. 



C. Rotation 

In this section we derive the response of the conden- 
sate phase and density to rotation. Without rotation, 
the wavefunction phase is constant for both components. 
After introducing rotation, solutions to the equations of 
motion can be found by inserting a quadrupolar oscilla- 
tion ansatz for the phase 

m 

4>, = -[a°H{p°)H{p°,) + a]H{-p°,)H{p])]xy, (8) 
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FIG. 1: (Color online) Examples of ID slices of the TF density 
profiles (dashed curves) and exact numerical solutions of the 
two component GPE (solid curves) for component 1 (black), 
component 2 (green [grey]), ei =62 = 0, cJi — UJ2 =0.1, 
h — I, mi — I, m2 — 1.5, gi — I, g2 — 0.5. (I) 5^12 = —0.5, 
(II) = 0.4, (III) gi2 = 0.65,(IV) gi2 = 0.8. The horizontal 
(vertical) axes are given in units of / (/~^), where / = gimi/h^. 

and Eq. (O into Eq. ([2]) and solving for the a's. For 
this gives: 

af - 2ap^ + ujj^ {a'j - Cj^) = 0, (9) 

which has up to 3 real solutions and is identical to the 
one component case ^]. For the a^'s we get two simul- 
taneous equations: 

-^2gjmf (a""/ - 2(x],^^ ^ uj [ex], - eyVi)^ = 0. (10) 

This yields up to nine real solutions for {aj,a2}. 

The a's are real constants representing the magnitude 
and orientation (positive or negative) of quadrupolar os- 
cillation in the two components. Specifically, the singular 
and overlapping regions of component 1 (component 2) 
undergo quadrupolar oscillations of different magnitudes, 
labeled and (al and 0^2) respectively. The density 
profile of each component is heavily influenced by this 
oscillation, in particular its ellipticity. This dependence 
has the interesting property that when a is negative, the 
BEC is deformed oppositely to the elliptical deformation 
caused by the trap. In other words the BEC forms an el- 
liptical profile that is rotated 90° to that of the confining 
elliptical potential. 

D. Discussion 

An example with 7 solutions to Eq. (p!Q|) for {af , 0^2} is 
shown in Fig.O Each branch represents a different static 
solution for a rotating TCBEC, each with quadrupolar 
oscillations of different magnitudes and orientations. Al- 
though these static solutions display a very complicated 
dependence on the many free parameters in the system. 
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lar oscillations in opposite directions (i.e. > and 
< 0). 

The BEC cannot follow these static solutions for all 
parameter values; eventually, one of two possibilities oc- 
curs. In the first case, the static solution which the BEC 
is following can cease to exist {catastrophic instability) 
[23 |. This leads to a massive disruption in the BEC's 
density profile, followed by the onset of turbulence and 
vortex nucleation. The second case occurs when the BEC 
is still in a state described by a static solution, but the 
solution itself is not stable. In this case, the BEC can 
either become turbulent {ripple instability) which leads 
to the formation of vortices, or the COM of the BEC 
becomes unstable {centre of mass instability). 



FIG. 2: (Color online) The solutions for the a's with H — 
111 — 112 — l,cJi = 0.2, a;2 = 0.1, ei — €2 — .l,mi = l,m2 = 
1.5,5'! = 1,^2 = 0.5,^12 = 0.1. Shown are solutions to ai 
(thick green [thick grey]), 0^2 (red dots [grey dots]), a\ (black 
dashed) , a2 (black line) . The branch numbering sits above the 
corresponding branch, it allows matching of each ai solutions 
to the corresponding solution. Q and a are in units of 

iii/h. 



the situation is greatly simplified by the fact that the 
BEC will prefer the branch with the lowest energy. 

The energy of a given solution increases as the mag- 
nitude of a increases. Also, for the solutions given by 
Eq. ([5]), the density of component i in the overlapping 
and singular regions does not join continuously at the 
boundary unless = This adds a large amount of 
kinetic energy at the interface of the two regions as the 
density profile must vary rapidly to join the two regions. 
The solution with the smallest energy can therefore be 
determined by removing branches that do not corre- 
spond closely to one of the branches. One then selects 
the remaining branch which has the smallest magnitude. 
From the figure, it is evident that for < 0.06 branch 
1 has the lowest energy, for 0.06 < Vt < 0.16 branch 2 
has the lowest energy, and for Q > 0.16 branch 4 has the 
lowest energy. 

The behaviour of the BEC can then be described as fol- 
lows. Initially, after condensation and before introducing 
rotation, the BEC will be in a state that corresponds 
to the stable static solution which has the lowest energy 
for the given set of experimental parameters. The BEC 
can then be transferred to a new state by adiabatically 
ramping these parameters. As this is done, the BEC 
will move along one of the static solutions. In the case 
of Fig. O each of branches 1, 2, and 4 can be accessed 
by ramping the parameters in different ways. Branch 1 
can be accessed by keeping e constant while ramping Q 
from to some final value or by keeping Q < 0.06 fixed 
and ramping e. Branch 2 (4) can be accessed by fixing 
0.06 < < 0.16 (O > 0.16) and ramping e from to some 
final value. Branch 2 is also of particular interest because 
it leads to each BEC component undergoing quadrupo- 



II. STABILITY 

A. Equations 

Stability can be analysed by linearising Eqs. (|2l3p 
about the critical points. We consider infinitesimal per- 
turbations pj pOj-\-Spj and (j)Oj -\-5(j)j ^ with {poi, 
(j)Oi} being a set of static solutions to Eqs. (|2|3p . In the 
overlapping region where both condensates coexist, one 
obtains 
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and in the singular region one obtains 
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Here vj = ^ V^Oj — ^ x i* is the wave function velocity in 
the rotating frame at position r. poj is given by Eq. ([5]) in 
the overlapping region and Eq. (|6]) in the singular region. 

As in the one component case [20j, we find that the 
eigenfunctions of the collective mode equations [Eqs. (pTI - 
[T3|) ] are polynomials of the form 5pj = ^pq^ fSjpqrX^y^z^ ^ 
^^3 = T.pqr^jpqrX^y'^z'', whcrc pjpqr and jjpqr are con- 
stants. The BEC is unstable when one of the eigenvalues 
has a positive real part, meaning that small perturbations 
about the static solutions grow exponentially. 

As well as the above method, extra information can be 
obtained by investigating the stability of the overlapping 
region of each component separately. This is done by sub- 
stituting Eq. ^ [instead of Eq. into Eq. By 
doing so, one can see in which component and in which 
region an instability originates. This, however, neglects 
interspecies cross terms; it is equivalent to treating the 
interaction between components as a static potential. In- 
terestingly, we find that this second method gives more 
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accurate results than if cross terms are considered. We 
take this to be an indication that the cohective mode 
cross interactions are not weh described within the TFA 
framework. 



B. Different types of instabilities 

Calculating the eigenvalues to Eq. (pT|) gives the re- 
gions of instabilitiy in the parameter space. Valuable in- 
formation regarding the nature of the instability can also 
be gained by observing the corresponding eigenvectors. 
This allows the instabilities of a rotating two-component 
BEC to be divided into four different types based on the 
very different effects they have on the overall structure 
of the BEC. Simulations showing the effects of these in- 
stabilties in detail are given in the next section. 



overlapping region of the condensate where both com- 
ponents are interacting. These intra- species COM insta- 
bilities are due purely to interactions of the superfluid 
components and lead to instability in the COM of each 
component separately, but the total COM of the conden- 
sate remains stable. They result in interesting dynamics 
which are described in section Hill 

3. Ripple Instability 

Perturbations of quadratic or higher in the position 
coordinates represent ripples through the phase and den- 
sity profile of the BEC. If these perturbations are un- 
stable they directly disrupt the smooth quadratic profile. 
These ripple instabilities lead to turbulence and vortex 
nucleation. 



1. Classical COM instability 

In its original static configuration, the density of BEC 
component j is a quadratic of the form (for simplicity 
of the argument we ignore trap anisotropy and the dis- 
tinction between the overlapping and singular region) 
Pj ~ 1^3 ~ ^'j ~ y] ~ ^'j' Consider the case where an 
eigenvalue A corresponds to an eigenvector first order 
in a position coordinate, e.g. one of the form Spj = 
Ax Xj . This perturbation will have the effect of displac- 
ing the centre of mass of the BEC along the x-axis, i.e. 
p Sp = {xj + .5Aa,)^ + ^1 + + fij. If A has a posi- 
tive real part then this perturbation will grow in magni- 
tude, displacing the centre of mass of the BEC even fur- 
ther and hence causing a COM instability. We find that 
component 1 always has such an instability in the range 
cjiv'l — ei < Q < cc^iV^l + ei regardless of the other pa- 
rameters. This is in fact the same instability as experi- 
enced by a classical point particle in a rotating harmonic 
trap. It is caused by the rotation frequency coupling 
to the oscillation frequency and is also experienced by a 
one-component BEC [19]. In this case, it either causes 
the BEC to oscillate as a whole about the trap centre, 
or it can drive the BEC out of the trap ^5[. In a two- 
component BEC, it can still occur provided there is little 
interaction between the components. This is discussed in 
more detail in Sec. Ill Ci 



2. Intra-species COM instability 

The above classical COM instability affects each com- 
ponent separately and independently of the interspecies 
interactions. We find another class of COM instability 
that occurs due to the interaction of the two components 
and can have a profoundly different effect on the BEC. 
They are again predicted by the instability of a pertur- 
bation first order in a position coordinate, but are dif- 
ferentiated from the first type in that they appear in the 



4. Catastrophic Instability 

In this case, the static solution that the BEC was fol- 
lowing during adiabatic ramping ceases to exist. Pertur- 
bations of all orders are unstable and the BEC is torn 
apart in a spectacular fashion [24]. After the initial on- 
set, the BEC becomes turbulent, and, as in the ripple 
instability, vortices nucleate. 



C. Results 

Using Eq. (pT]) we present two examples (Figs. [3] and 
m of the instabilities predicted by the TFA and how they 
appear in different regions of the parameter space. We 
have used the second method described in Sec. [Ill On 
the same figures are plotted the point where the BEC 
becomes unstable in GPE simulations conducted for the 
same set of parameters. In general the solutions and their 
stability can be evaluated for any ramping procedure, the 
only difference being that different ramping paths can 
accesses different a branches. In our GPE simulations 
we fix Q and then ramp ei = 62 from 0. This allows 
us to investigate regions of the phase diagram that, for 
example, would not be accessible by ramping Q for fixed 
e. The point of instability is then determined from the 
simulations using the method presented in [25|. 

The results of the simulations show how the TFA re- 
sults can be used to predict BEC behaviour. The meth- 
ods derived above successfully pick up the different re- 
gions in the parameter space where different instabilities 
occur. Sometimes the results of the GPE simulations are 
exactly predicted in the TFA, at other times the phe- 
nomenology is as predicted but its location numerically 
shifted in ft by as much as 15% . As discussed above, 
the TFA is valid within the bulk of the condensate. The 
shifting is an indication that for the chosen parameters, 
the boundary of the BEC is affecting the particular in- 
stability that has been shifted. 
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For the case of repulsive interspecies interaction {gi2 > 
0), we find that the stability of the overlapping region of 
component 1 is highly dependent on its connection to the 
singular region [see, for example, the non rotating case 
Fig.[T](III)]. The results using the TEA on the overlapping 
region of component 1 are only accurate when gi2 < 
and should not be used for gi2 > 0. 

Fig. [3] shows the instabilities found on branch 2 of Fig. 
[21 A BEC on this branch can, depending on the pa- 
rameters, experiences either a intra-species COM, ripple, 
catastrophic, or classical COM instability. The actual po- 
sition of the catastrophic instability is shifted to higher Q 
from the TFA prediction by ~ 10% and the intra-species 
COM instability is shifted to higher Qhy 5%. The green 
dashed line indicates the region where the classical COM 
instability would be for component 2 were the other com- 
ponent not present. 

Fig. m focuses on the COM instability for a different 
system. This system has attractive interspecies inter- 
actions and equal masses but different trapping frequen- 
cies, which could correspond, for example, to two conden- 
sates of the same atomic species but different hyperfine 
states. The intra-species COM instability and ripple in- 
stability in component 2 are exactly as predicted. The 
catastrophic instability is shifted to lower Q from the pre- 
dictions by 10%. The ripple instability in the singular 
region of component 1 is similar to, but much less influen- 
tial than, the prediction. However this is not surprising 
because for attractive interspecies interactions the two 
components are pulled tightly together and the singular 
region is almost non-existent [see Fig. [T](I)]. 

The TFA predictions for both examples show that the 
interspecies interactions should halt the onset of the clas- 
sical COM instabilities when both species are performing 
stable quadrupolar oscillation. Interestingly, the intra- 
species COM instability pulls the two components apart 
leaving them once again susceptible to the classical COM 
instability. If this has occurred before reaching, or whilst 
within the classical COM instability, the component in 
question will exit the trap as though the other component 
were not present. Likewise, the catastrophic instability 
causes both components to be wrenched apart. If one 
of the components is within a classical COM instability 
during the break down, it is free of the other component 
for long enough that it will exit the trap. This effect is 
clearly seen on Figs [3] and [H 



III. SIMULATIONS 

To compliment these Thomas-Fermi results we have 
adapted 2D Truncated Wigner simulations [s?] to the 
two-component case, and investigated the instabilities 
predicted. Interesting results are attained for each dif- 
ferent type of instability, and these can be understood in 
terms of the TFA results. Also, the numerical methods 
allow us to extend the investigation to the immiscible 
phase where the TFA is no longer valid. 
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FIG. 3: (Color online) The instability regions within the pa- 
rameter space for a system with h — iii —112 = l,<^i — 
0.2, mi — l,gi — 1,CJ2 — 0.1, 7712 — 1.5,^2 — 0.5, 5^12 — .1, 
catastrophic instability (Region I - black), intra-species COM 
instability (Region II - red [dark grey]), ripple instability in 
component 2 (Region III - orange [light grey]), ripple insta- 
bility in component 1 (Region IV - blue [dark grey]), classical 
COM instability in component 2 (green [light grey] dashed 
line). Marks show the results for instabilities found using 
GPE simulations for catastrophic instabilities (circles), intra- 
species COM instabilities (squares), classical COM instabili- 
ties (triangles), ripple instabilities (crosses). Q is in units of 

iii/h. 
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FIG. 4: (Color online) same as Fig.[3]but with = /ii = /i2 = 
l,a;i = 0.1, mi = l^gi — 0.5, a;2 = 0.15, m2 = 1,5'2 = l,5'i2 = 
-.3. 



As well as the results shown in Figs. [3] and [4] we conduct 
simulations for parameters that correspond to an ^^Rb- 
i^^Cs system [ss] and an ^'^Rb-^^Rb system [39]. These 
systems are of particular interest because the scattering 
lengths of ^^^Cs and ^^Rb can be tuned via a Feshbach 
resonance [SO*, ^4^. This means that both the immisci- 
ble and miscible phases can be accessed experimentally. 
The ^^Rb-^^Rb is also important because it satisfies the 
requirements for the creation of vortex sheets. The ^^Rb- 
^^^Cs system is used to simulate the creation of the other 
rotating states. For the ^^Rb-^^^Cs system, ai2 is un- 
known. We, therefore, give examples of simulations con- 
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ducted for a number of different choices of ai2. 



1. The truncated Wigner method 

The truncated Wigner method simulates quantum vac- 
uum fluctuations by adding appropriate classical random 
fluctuations to the coherent field of the EEC's initial 
state. In this system the fluctuations serve two pur- 
poses. Firstly, the fluctuations provide a seed of noise to 
break the BEC symmetry when the quadrupolar oscilla- 
tion becomes unstable. Secondly, they enable incoherent 
scattering processes to occur, by which condensate atoms 
are scattered into a thermal cloud. This method has been 
used to describe, for example, the formation of scattering 
halos in condensate collisions ^Hj, [42[, and the suppres- 
sion of Cherenkov radiation [43|. Hence the turbulent 
BEC can relax into a rotating eigenstate, such a vortex 
lattice, in contrast to the bare GPE which conserves en- 
ergy and atom number. 

In practice, the fluctuations are included as follows. 
The initial wavefunction for each species is obtained 
by solving the time-independent two-component GPE. 
These two wavefunctions are then expanded over a plane- 
wave basis, with a maximum cutoff wave vector to pre- 
vent Fourier aliasing. Quantum fluctuations are intro- 
duced into each wavefunction separately by adding ran- 
dom complex noise to each plane- wave mode. The ampli- 
tude of the quantum fluctuations has a Gaussian distri- 
bution, with an average value of half a particle [4lJ . For 
a thorough description of the method, see Refs. [33, [42[. 
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FIG. 5: (Color online) Plots of the density for BECs undergo- 
ing intra-species COM instabilities for an ^^Rb (component 1 
- blue [panels Bl B2, right half of Al and bottom half of A2]) - 
^^^Cs (component 2 - red [panels CI C2, left half of Al and top 
half of A2]) system with ai = a2 = 5.4 nm, cji = a;2 = 3.15 
Hz. Attractive case {gi2 < 0): ai2 = —0.65 nm, Q = 3.66 
Hz. (A2) is the density shortly after (Al) showing the two 
components rotating around each other. The black dot is the 
COM of the system. Repulsive case {gi2 > 0): ai2 = 3.9 nm, 
Q = 3.66 Hz. (Bl, CI) is the density shortly after (B2,C2) 
showing component 2 (C1,C2) bouncing of the shell of the 
enclosing component 1 (B1,B2). The horizontal bars denote 
40 /im. 



settle down into a state that has component 2 still sitting 
within component 1, but with both components signifi- 
cantly diffused into one another. During these instabili- 
ties, no vortices are formed. 



A. Results 

1. Intra-species COM instability 

The TEA results show that as gi2 is increased or de- 
creased from 0, the classical COM instability of com- 
ponent 2 begins to shift and becomes an intra-species 
COM instability. For gi2 < the simulations show that 
when this instability is reached, the COM of the entire 
BEC is stable; however, the COM of the individual com- 
ponents becomes unstable and the two components sep- 
arate. The instability is caused purely by interspecies 
interactions, so as the components separate and interac- 
tions decrease, the trap pushes the two components back 
together again. The effect is that both components be- 
gin oscillating; eventually, the BEC settles down into a 
stable state with both components orbiting around one 
another [Fig. [5] (A)]. This state is static in the rotating 
frame and implies a breaking of the 180° symmetry of 
the rotating frame Hamiltonian. 

For gi2 > component 2 sits within component 1. 
When the intra-species COM instability in component 
2 is reached, it is still trapped within component 1; it 
begins bouncing off the surrounding shell [Fig. [5] (B,C)]. 
Eventually, the two components disrupt each other and 



2. Classical COM instability 

The classical COM instability can affect a TCBEC if 
both components experience it at the same time (i.e. if 
LJi ~ 6^2). However, if only one component is in a regime 
of classical COM instability then the TFA results pre- 
dict, and simulations confirm, that this instability is sup- 
pressed during stable rotational motion. None-the-less, 
as already discussed in Sec. Ill Cl each component can still 
experience its classical COM instability independently of 
the other component provided it is freed from the other 
component first via a catastrophic or intra-species COM 
instability. 

Another case where a classical COM instability can 
occur is near or in the immiscible phase where both com- 
ponents have a large singular region. In this case, com- 
ponent 2 forms a tightly packed ball within the shell of 
component 1. When component 2 reaches its classical 
COM instability it builds up enough momentum to break 
through component 1 while component 1 remains stable. 
Component 2 then either exits the trap or begins to oscil- 
late within the trapping potential. The result is similar to 
a wrecking ball as the highly dense component 2 smashes 
through the dilute component 1. 
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FIG. 6: (Color online) Results of 2D Truncated Wigner simulations ending in vortex nucleation. Panels D,E,F,and G show 
an ^^Rb (component 1) - "^^^Cs (component 2) system with ai = a2 = 5.4 nm, u;i = a;2 = 3.15 Hz. (D,E) ai2 = 7 nm, (F,G) 
ai2 = 3.3 nm. Panels H and I show an ^^Rb (component 1) - ^^Rb (component 2) system with ai = 5.24 nm, a2 = 11.27 
nm, ai2 = 11.3 nm, cji = a;2 = 3.15 Hz. The horizontal bars denote 40 /xm. The number of atoms in each simulation is such 
that the peak density of ^^Rb is 4 x 10^^ m~^ and both components have equal norms. D1,F1, and HI (E1,G1, and II) is the 
density and D2,F2, and H2 (E2,G2, and 12) the phase of component 1 (2). The simulations are conducted by ramping e while 
fixing Q at 1.75 Hz 



3. Catastrophic and ripple instabilities 

In the miscible phase, the simulations show that when 
a ripple or catastrophic instability is reached the BEG 
becomes turbulent and vortices enter. They then form 
an interlaced vortex lattice [Fig. [6] (F,G)]. This state is 
of the same type as has been seen to form in experiments 
when a one-component BEG with a vortex lattice already 
present was split into two hyperfine components [i3| . 

We also investigated the behaviour of the rotational 
instabilities near to and within the immiscible phase. 
When gi2 is close to but less than ^/gigi^ component 
1 forms an almost hollow shell around component 2 even 
though the immiscible phase has not been reached. As 
expected from the discussion at the end of Sec. II Bt once 
the condensate is stirred and turbulence reached, the two 
components mix together and a vortex lattice is formed. 

Once within the immiscible phase however, there is 
a discontinuous change in the behaviour of a rotated 
BEG, and, importantly, the two components do not mix. 
We find that once turbulence is induced, two different 
possible configurations with non trivial phase topology 
form. The first is a vortex sheet which has been shown 
to be preferable thermodynamically in regions of param- 
eter space where gi ^ ^2, ^1 ^ ^2, and gi2 > ^Jglg2 
|l5| [Fig. [6] (H, I)]. In this configuration the two com- 
ponents form separate domains of high density. We find 
that the parameters required for the dynamical forma- 
tion of a vortex sheet matches that of the thermodynamic 
analysis provided an instability has been reached; with- 
out an instability they do not form. These conditions 
for vortex sheet formation also coincide with the param- 
eters required for a symmetry broken initial state in the 
absence of rotation [Sq . 

The second state we find is a giant vortex which has 
also been predicted thermodynamically [IT,, il8t] . We find 



that this state forms for any set of parameters that fall 
within the immiscible phase but do not match the re- 
quirements for the formation of vortex sheets. As with 
the vortex sheets, a ripple or catastrophic instability 
must first be induced in order for this state to form dy- 
namically. The simulations show that once instability 
is reached multiple vortices push their way through the 
outer shell of component one and into the low density 
region in the middle of the trap where they congregate 
forming a giant vortex [Fig.[6](D,E)]. Normally, two over- 
lapping vortices are predicted to be thermodynamically 
unstable [44]. However, in this system the large density of 
component 2 attracts them to the centre while the outer 
shell of Kb holds them in. These results show how in 
practice such a state could be created and that this state 
is a product of the immiscibility of the two components. 

In summary, the above results show that there is a 
direct connection between the state of the non-rotating 
TGBEG and the state the rotating TGBEG will settle 
down into after instability has been reached. That is, 
a miscible phase leads to vortex lattice after rotating, a 
symmetry broken state leads to vortex sheets, and an im- 
miscible phase (other than the symmetry broken phase) 
leads to a giant vortex. 



IV. CONCLUSION 

These results illustrate that BEG mixtures produce a 
rich variety of dynamical regimes that may be accessed 
by tuning experimental parameters. We have examined 
the applicability of the TFA to a TGBEG and extended it 
to the rotating case allowing us to find symmetry break- 
ing GOM oscillations that are induced by interspecies 
interactions. In addition, the results give conditions un- 
der which interlaced vortex lattices, giant vortices and 
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vortex sheets will spontaneously form. The method al- the EPSRC (RGS). 
lows all these phenomena to be understood and classi- 
fied through particular instabilities occurring in different 
parts of the total condensate. 
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